home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 376-400 / disk_386 / xlispstat / src2.lzh / XLisp-Stat / lowess.c < prev    next >
C/C++ Source or Header  |  1990-10-04  |  5KB  |  169 lines
  1. /*Translated from RATFOR lowess code of W. S. Cleveland as obtained from NETLIB */
  2. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  3. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  4. /* You may give out copies of this software; for conditions see the    */
  5. /* file COPYING included with this distribution.                       */
  6.  
  7. #include <stdlib.h>
  8. #include "xlisp.h"
  9. #include "osdef.h"
  10. #ifdef ANSI
  11. #include "xlsproto.h"
  12. #else
  13. #include "xlsfun.h"
  14. #endif ANSI
  15.  
  16. #ifdef ANSI
  17. double pow2(double),pow3(double),fmax(double,double);
  18. void lowest(double *,double *,int,double,double *,int,int,double *,int,double *,
  19.      int *),sort(double *,int);
  20. int compar(double *,double *);
  21. #else
  22. double pow2(),pow3(),fmax();
  23. void lowest(),sort();
  24. int compar();
  25. #endif ANSI
  26.  
  27. #define FALSE 0
  28. #define TRUE 1
  29.  
  30. static double pow2(x) double x; { return(x * x); }
  31. static double pow3(x) double x; { return(x * x * x); }
  32. static double fmax(x,y) double x, y; { return (x > y ? x : y); }
  33.  
  34. int lowess(x, y, n, f, nsteps, delta, ys, rw, res)
  35.      double /* *x, *y,*/ f, delta/*, *ys, *rw, *res*/;/* changed JKL */
  36.      RVector x, y, ys, rw, res;
  37.      int n, nsteps;
  38. {
  39.   int iter, ns, ok, nleft, nright, i, j, last, m1, m2;
  40.   double d1, d2, denom, alpha, cut, cmad, c9, c1, r;
  41.   
  42.   if (n < 2) { ys[0] = y[0]; return(1); }
  43.   ns = max(min((int) (f * n), n), 2);  /* at least two, at most n points */
  44.   for(iter = 1; iter <= nsteps + 1; iter++){      /* robustness iterations */
  45.     nleft = 0; nright = ns - 1;
  46.     last = -1;        /* index of prev estimated point */
  47.     i = 0;   /* index of current point */
  48.     do {
  49.       while(nright < n - 1){
  50.     /* move nleft, nright to right if radius decreases */
  51.     d1 = x[i] - x[nleft];
  52.     d2 = x[nright + 1] - x[i];
  53.     /* if d1 <= d2 with x[nright+1] == x[nright], lowest fixes */
  54.     if (d1 <= d2) break;
  55.     /* radius will not decrease by move right */
  56.     nleft++;
  57.     nright++;
  58.       }
  59.       lowest(x, y, n, x[i], &ys[i], nleft, nright, res, (iter > 1), rw, &ok);
  60.       /* fitted value at x[i] */
  61.       if (! ok) ys[i] = y[i];
  62.       /* all weights zero - copy over value (all rw==0) */
  63.       if (last < i - 1) { /* skipped points -- interpolate */
  64.     denom = x[i] - x[last];    /* non-zero - proof? */
  65.     for(j = last + 1; j < i; j = j + 1){
  66.       alpha = (x[j] - x[last]) / denom;
  67.       ys[j] = alpha * ys[i] + (1.0 - alpha) * ys[last];
  68.     }
  69.       }
  70.       last = i;        /* last point actually estimated */
  71.       cut = x[last] + delta;     /* x coord of close points */
  72.       for(i=last + 1; i < n; i++) {     /* find close points */
  73.     if (x[i] > cut) break;     /* i one beyond last pt within cut */
  74.     if(x[i] == x[last]) {      /* exact match in x */
  75.       ys[i] = ys[last];
  76.       last = i;
  77.     }
  78.       }
  79.       i = max(last + 1,i - 1);
  80.       /* back 1 point so interpolation within delta, but always go forward */
  81.     } while(last < n - 1);
  82.     for (i = 0; i < n; i++)      /* residuals */
  83.       res[i] = y[i] - ys[i];
  84.     if (iter > nsteps) break; /* compute robustness weights except last time */
  85.     for (i = 0; i < n; i++) 
  86.       rw[i] = fabs(res[i]);
  87.     sort(rw,n);
  88.     m1 = 1 + n / 2; m2 = n - m1 + 1;
  89.     cmad = 3.0 * (rw[m1] + rw[m2]);      /* 6 median abs resid */
  90.     c9 = .999 * cmad; c1 = .001 * cmad;
  91.     for (i = 0; i < n; i++) {
  92.       r = fabs(res[i]);
  93.       if(r <= c1) rw[i] = 1.0;      /* near 0, avoid underflow */
  94.       else if(r > c9) rw[i] = 0.0;  /* near 1, avoid underflow */
  95.       else rw[i] = pow2(1.0 - pow2(r / cmad));
  96.     }
  97.   }
  98.   return(0);
  99. }
  100.  
  101.  
  102. static void lowest(x, y, n, xs, ys, nleft, nright, w, userw, rw, ok)
  103.      double *x, *y, *w, *rw, xs, *ys;
  104.      int n, nleft, nright, userw, *ok;
  105. {
  106.   double range, h, h1, h9, a, b, c, r;
  107.   int j, nrt;
  108.  
  109.   range = x[n - 1] - x[0];
  110.   h = fmax(xs - x[nleft], x[nright] - xs);
  111.   h9 = .999 * h;
  112.   h1 = .001 * h;
  113.  
  114.   /* compute weights (pick up all ties on right) */
  115.   a = 0.0;        /* sum of weights */
  116.   for(j = nleft; j < n; j++) {
  117.     w[j]=0.0;
  118.     r = fabs(x[j] - xs);
  119.     if (r <= h9) {    /* small enough for non-zero weight */
  120.       if (r > h1) w[j] = pow3(1.0-pow3(r/h));
  121.       else w[j] = 1.0;
  122.       if (userw) w[j] = rw[j] * w[j];
  123.       a += w[j];
  124.     }
  125.     else if (x[j] > xs) break;  /* get out at first zero wt on right */
  126.   }
  127.   nrt = j - 1;  /* rightmost pt (may be greater than nright because of ties) */
  128.   if (a <= 0.0) *ok = FALSE;
  129.   else { /* weighted least squares */
  130.     *ok = TRUE;
  131.  
  132.     /* make sum of w[j] == 1 */
  133.     for (j = nleft; j <= nrt; j++) w[j] = w[j] / a;
  134.  
  135.     if (h > 0.0) {     /* use linear fit */
  136.  
  137.       /* find weighted center of x values */
  138.       for (j = nleft, a = 0.0; j <= nrt; j++) a += w[j] * x[j];
  139.  
  140.       b = xs - a;
  141.       for (j = nleft, c = 0.0; j <= nrt; j++) 
  142.     c += w[j] * (x[j] - a) * (x[j] - a);
  143.  
  144.       if(sqrt(c) > .001 * range) {
  145.     /* points are spread out enough to compute slope */
  146.     b = b/c;
  147.     for (j = nleft; j <= nrt; j++) 
  148.       w[j] = w[j] * (1.0 + b*(x[j] - a));
  149.       }
  150.     }
  151.     for (j = nleft, *ys = 0.0; j <= nrt; j++) *ys += w[j] * y[j];
  152.   }
  153. }
  154.  
  155. static int compar(a, b)
  156.      double *a,*b;
  157. {
  158.   if (*a < *b) return(-1);
  159.   else if (*a > *b) return(1);
  160.   else return(0);
  161. }
  162.  
  163. static void sort(x, n)
  164.      double *x;
  165.      int n;
  166. {
  167.   qsort((char *)x, n, sizeof(double), compar); /* cast added JKL */
  168. }
  169.